Load packages

library(amerifluxr) # Package for querying, downloading, and handling AmeriFlux data and metadata.
library(REddyProc)  # Package for processing Eddy Covariance flux data (e.g., gap-filling, partitioning)
library(lubridate)  # Provides tools to work with date and time (e.g., parsing, formatting).
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(tidyverse)  # A collection of packages for data manipulation, visualization, and more.
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr   1.1.4     ✔ readr   2.1.5
## ✔ forcats 1.0.0     ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1     ✔ tibble  3.2.1
## ✔ purrr   1.0.2     ✔ tidyr   1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(bigleaf)  # Provides functions to calculate leaf area index and other plant trait estimations.
library(data.table)  # Offers fast and memory-efficient manipulation of large datasets.
## 
## Attaching package: 'data.table'
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
## 
## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
library(ggplot2)

01: Download AmeriFlux BASE publish

# check out data coverage for different sites
full.list = amf_data_coverage(data_product = "BASE-BADM", data_policy = "CCBY4.0")
head(full.list)
##   SITE_ID                                             URL
## 1  AR-Bal https://ameriflux.lbl.gov/sites/siteinfo/AR-Bal
## 2  AR-CCa https://ameriflux.lbl.gov/sites/siteinfo/AR-CCa
## 3  AR-CCg https://ameriflux.lbl.gov/sites/siteinfo/AR-CCg
## 4  AR-Cel https://ameriflux.lbl.gov/sites/siteinfo/AR-Cel
## 5  AR-Rom https://ameriflux.lbl.gov/sites/siteinfo/AR-Rom
## 6  AR-TF1 https://ameriflux.lbl.gov/sites/siteinfo/AR-TF1
##                                          publish_years
## 1                                           2012, 2013
## 2 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020
## 3             2018, 2019, 2020, 2021, 2022, 2023, 2024
## 4                                                     
## 5                                                     
## 6                                     2016, 2017, 2018
print(paste("number of sites included in BASE publish is", nrow(full.list)))
## [1] "number of sites included in BASE publish is 759"
# download single site flux/met data
site_name = "US-Syv"

# please keep the codes below commented if you do not have an account 
# floc<- amf_download_base(
#     user_id = "my_user", # change to your user name
#     user_email = "my_email@gmail.com", # change to your email
#     site_id = site_name,
#     data_product = "BASE-BADM",
#     data_policy = "CCBY4.0",
#     agree_policy = TRUE,
#     intended_use = "other",
#     intended_use_text = "Flux-LSM_workshop",
#     verbose = TRUE,
#     out_dir = tmp_dir
#   )

# If you do not have an account, you can directly start from here
tmp_dir <- tempdir()
setwd("~/") # binder wd
setwd("C:/Users/yl763/Documents/GitHub/FCC_workshop_flux_test") # Yujie's local wd, will be removed later!
floc =  "AMF_US-Syv_BASE-BADM_29-5.zip"
base <- amf_read_base(file = floc,
                      unzip = TRUE,
                      parse_timestamp = TRUE) 

# check data coverage for meteorology
# visualizes the BASE data availability for selected AmeriFlux sites, variables, and years.
##PPFD_IN   Photosynthetic photon flux density, incoming    µmolPhoton m-2 s-1
## TA   Air temperature deg C
## RH   Relative humidity, range 0-100  %
## VPD  Vapor Pressure Deficit  hPa
# more variable names: https://ameriflux.lbl.gov/data/aboutdata/data-variables/#base
amf_plot_datayear(site_set = site_name, 
                  var_set = c("PPFD_IN", "TA", "RH", "VPD"),
                  # year_set = c(2004:2023),
                  nonfilled_only = TRUE)
# check data coverage for flux data
amf_plot_datayear(site_set = site_name, 
                  var_set = c("FC", "H", "LE"),
                  # year_set = c(2004:2023),
                  nonfilled_only = TRUE)
# You can uncomment the line below to check data availability by year for all variables
# amf_plot_datayear(site_set = site_name, nonfilled_only = FALSE)

02: Organise input data

# subset data for two years
base = base[base$YEAR %in% c(2020,2021), ]

# Recreate other time variables using TIMESTAMP (This is just because I sometimes get error messages in the next step "Intialize EProc" when using the default time variables.)
recreate_time_vars <- function(df) {
  df %>%
    mutate(
      TIMESTAMP = substr(TIMESTAMP_START, 1, 12),
      year = substr(TIMESTAMP_START, 1, 4),
      month = substr(TIMESTAMP_START, 5, 6),
      day = substr(TIMESTAMP_START, 7, 8),
      hour = as.numeric(substr(TIMESTAMP_START, 9, 10)) + c(0.5, 1),
      date = as.Date(paste(year, month, day, sep = "-")),
      doy = yday(date)
    )
}
base <- recreate_time_vars(base)
# format the data to be used as input for REddyProc
base_df = data.frame(
  'Year' = base$YEAR,
  'Hour' = base$hour,
  'Date' = base$date,
  'Month' = base$month,
  'DoY' = base$DOY,
  'NEE' = base$FC,
  'LE' = base$LE,
  'H' = base$LE,
  'Rg' = ifelse(base$PPFD_IN_PI_F_1_1_1 < 0, 0, PPFD.to.Rg(base$PPFD_IN_PI_F_1_1_1)),
  'Tair' = base$TA_1_1_1,
  'Tsoil' = base$TS_1_1_1,
  'RH' = ifelse(base$RH_1_1_1 > 100, 100, base$RH_1_1_1),
  'VPD' = base$VPD_PI_F_1_1_1,
  'Ustar' = base$USTAR_1_1_1,
  'TIMESTAMP_START' = as.character(base$TIMESTAMP_START),
  'TIMESTAMP_END' = as.character(base$TIMESTAMP_END),
  'PPFD' = base$PPFD_IN_PI_F_1_1_1
)

head(base_df)
##   Year Hour       Date Month DoY       NEE        LE         H Rg      Tair
## 1 2020  0.5 2020-01-01    01   1 0.6667798 -4.066935 -4.066935  0 -12.74600
## 2 2020  1.0 2020-01-01    01   1        NA        NA        NA  0 -13.31267
## 3 2020  1.5 2020-01-01    01   1        NA -2.449585 -2.449585  0 -13.09900
## 4 2020  2.0 2020-01-01    01   1        NA        NA        NA  0 -12.97467
## 5 2020  2.5 2020-01-01    01   1 0.5290407 -2.824548 -2.824548  0 -13.11433
## 6 2020  3.0 2020-01-01    01   1        NA        NA        NA  0 -13.14767
##       Tsoil       RH       VPD     Ustar TIMESTAMP_START TIMESTAMP_END PPFD
## 1 0.6363515 87.22335 0.2930273 0.2653910    202001010000  202001010030    0
## 2 0.6407161 88.82667 0.2447449 0.2953020    202001010030  202001010100    0
## 3 0.6424356 89.99333 0.2230272 0.3732584    202001010100  202001010130    0
## 4 0.6455239 90.04335 0.2241620 0.6959635    202001010130  202001010200    0
## 5 0.6453573 89.80999 0.2268308 0.8161998    202001010200  202001010230    0
## 6 0.6477322 89.56665 0.2316199        NA    202001010230  202001010300    0

03: Intialize EProc

?filterLongRuns
## starting httpd help server ... done
#filterLongRuns : Longer runs, i.e. sequences of numerically identical values, in a series of measurements hint to problems during a noisy measurement, e.g. by sensor malfunction due to freezing. This function, replaces such values in such runs with NA to indicate missing values.
EddyData <- filterLongRuns(base_df, "NEE")
EddyData$Year <- as.numeric(EddyData$Year)
EddyData$Hour <- as.numeric(EddyData$Hour)
EddyData$DoY <- as.numeric(EddyData$DoY)
EddyDataWithPosix <- fConvertTimeToPosix(EddyData, 'YDH', Year = 'Year', Day = 'DoY', Hour  = 'Hour') 
## Converted time format 'YDH' to POSIX with column name 'DateTime'.
EProc <- sEddyProc$new(site_name, EddyDataWithPosix, c('NEE','Rg','Tair','VPD', 'Ustar', "H", "LE"))
## New sEddyProc class for site 'US-Syv'
# EProc an R object used in REddyProc with the attributes defined
class(EProc)
## [1] "sEddyProc"
## attr(,"package")
## [1] "REddyProc"

04: IQR filtering

# Gapfill meteorological variables using MDS
EProc$sMDSGapFill('Tair', FillAll = FALSE,  minNWarnRunLength = NA)
## Initialized variable 'Tair' with 2810 real gaps for gap filling.
## Limited MDS algorithm for gap filling of 'Tair' with LUT(Rg only) and MDC.
## Look up table with window size of 7 days with Rg
## ............................1568
## Mean diurnal course with window size of 0 days: .
## ............0
## Mean diurnal course with window size of 1 days: .
## ............0
## Mean diurnal course with window size of 2 days: .
## ............0
## Look up table with window size of 14 days with Rg
## ............667
## Look up table with window size of 21 days with Rg
## .....575
## Finished gap filling of 'Tair' in 1 seconds. Artificial gaps filled: 35088, real gaps filled: 2810, unfilled (long) gaps: 0.
EProc$sMDSGapFill('Rg', FillAll = FALSE,  minNWarnRunLength = NA)
## Initialized variable 'Rg' with 0 real gaps for gap filling.
## Restriced MDS algorithm for gap filling of 'Rg' with no meteo conditions and hence only MDC.
## Finished gap filling of 'Rg' in 0 seconds. Artificial gaps filled: 35088, real gaps filled: 0, unfilled (long) gaps: 0.
EProc$sMDSGapFill('VPD', FillAll = FALSE,  minNWarnRunLength = NA)
## Initialized variable 'VPD' with 4270 real gaps for gap filling.
## Limited MDS algorithm for gap filling of 'VPD' with LUT(Rg only) and MDC.
## Look up table with window size of 7 days with Rg
## ..........................................2422
## Mean diurnal course with window size of 0 days: .
## ..................0
## Mean diurnal course with window size of 1 days: .
## ..................0
## Mean diurnal course with window size of 2 days: .
## ..................0
## Look up table with window size of 14 days with Rg
## ..................701
## Look up table with window size of 21 days with Rg
## ...........673
## Look up table with window size of 28 days with Rg
## ....474
## Finished gap filling of 'VPD' in 1 seconds. Artificial gaps filled: 35088, real gaps filled: 4270, unfilled (long) gaps: 0.
# Use MDS to get the "best guess" of true flux
EProc$sMDSGapFill('NEE') 
## Initialized variable 'NEE' with 14679 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29511
## Look up table with window size of 14 days with Rg VPD Tair
## .......................................................582
## Look up table with window size of 7 days with Rg
## .................................................2982
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................770
## Look up table with window size of 21 days with Rg
## ...........685
## Look up table with window size of 28 days with Rg
## .....513
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 14679, unfilled (long) gaps: 0.
# Calculate residuals and identify outliers
residual <- EProc$sTEMP$NEE_orig - EProc$sTEMP$NEE_fall
IQR <- IQR(residual, na.rm = TRUE)
outlier <- ifelse(abs(residual) > (IQR * 6), 1, 0)
EddieC <- data.frame(
    sDateTime = EProc$sTEMP$sDateTime,
    NEE_orig = EProc$sTEMP$NEE_orig,
    Ustar = EProc$sDATA$Ustar,
    NEE_fall = EProc$sTEMP$NEE_fall,
    residual = residual,
    outlier = outlier
  )
  
# Rename columns
colnames(EddieC) <- c('sDateTime', 'NEE_orig', 'Ustar', 'NEE_fall', 'residual', 'outlier')
  
# Filter out outliers
EddieC$NEE_filt <- dplyr::if_else(EddieC$outlier > 0, NA_real_, EddieC$NEE_orig)
EddieC$year <- substr(EddieC$sDateTime, 1, 4)
EddieC$doy <- strftime(EddieC$sDateTime, format = "%j")
# Plot the initial outlier detection
EddieC %>%
    arrange(as.factor(outlier)) %>%
    ggplot(aes(y = NEE_orig, x = as.numeric(doy), color = as.factor(outlier))) +
    geom_point(shape = 20, alpha = 0.4) +
    theme_minimal() +
    labs(x = 'Day of year', y = 'NEE') +
    scale_color_manual(values = c('skyblue', 'red')) +
    facet_wrap(~year) + 
    ylim(c(-50, 50)) +
    ggtitle("intial outlier detection")
## Warning: Removed 14679 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Re-run the outlier test after initial filtering
EddieC$residual2 <- EddieC$NEE_filt - EddieC$NEE_fall
EddieC$IQR2 <- IQR(EddieC$residual2, na.rm = TRUE)
EddieC$outlier2 <- ifelse(abs(EddieC$residual2) > EddieC$IQR2 * 6, 1, 0)
EddieC$NEE_filt2 <- ifelse(EddieC$outlier2 == 0, EddieC$NEE_filt, NA)
  
# Plot the re-run outlier detection
EddieC %>%
    arrange(as.factor(outlier2)) %>%
    ggplot(aes(y = NEE_orig, x = as.numeric(doy), color = as.factor(outlier2))) +
    geom_point(shape = 20) +
    theme_minimal() +
    labs(x = 'Day of year', y = 'NEE') +
    scale_color_manual(values = c('skyblue', 'red')) +
    facet_wrap(~year) + 
    ylim(c(-50, 50)) +
    ggtitle("re-run outlier detection")
## Warning: Removed 14679 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Remove outliers from the main data
EProc$sDATA$NEE <- EddieC$NEE_filt2

05: u* filtering

What is u*?

  • u*: Friction velocity (m/s) is a measure of turbulent mixing in the lower atmosphere near the surface;
  • It can be thought of as a velocity scale, a representative value for a ‘turbulent velocity’;
  • u* is defined using the vertical flux of horizontal momentum;
  • (u*)^2 = overbar(u’w’)^2 + overbar(v’w’)^2;
  • CO2 transport is heavily influenced by turbulent mixing;
  • Unfavorable conditions could be detected by inspecting the relationship of nighttime NEE vs. u*.

Why is u* filtering needed?

  • To filter out period with low turbulence;
  • At low u* values, respiration is negatively biased;
  • At night when winds are calm, the atmosphere often becomes stable with low turbulence and low mixing locally;
  • This could lead to CO2 pooling near the surface or within canopies;
  • EC towers may underestimate respiration fluxes as CO2 is not being transported upward efficiently to be measured.

Different u* threshold treatment options:

The u* threshold is the minimum u* above which respiration reaches a plateau (see figure under /plots). - User-specific u* threshold
- Single (fixed) u* threshold
- Annually varying u* threshold
- Seasonally varying u* threshold - More details: REddyProc u* cases vignette

set.seed(2000)
# Here, we are using nSample = 10L for demonstration, please use nSample = 1000L for real research 
uStarTh <- EProc$sEstUstarThresholdDistribution(nSample = 10L, probs = c(0.05, 0.5, 0.95)) 
## Warning in .estimateUStarSeason(...): sEstUstarThreshold: too few finite
## records within season (n = 622) for 7 temperature classes. Need at least n =
## 700. Returning NA for this Season.
## Warning in .estimateUStarSeason(...): sEstUstarThreshold: too few finite
## records within season (n = 633) for 7 temperature classes. Need at least n =
## 700. Returning NA for this Season.
## 
## Estimated UStar distribution of:
##        uStar        5%       50%       95%
## 1 0.4541172 0.4188225 0.4564736 0.5286447 
## by using  10 bootstrap samples and controls:
##                        taClasses                    UstarClasses 
##                               7                              20 
##                           swThr            minRecordsWithinTemp 
##                              10                             100 
##          minRecordsWithinSeason            minRecordsWithinYear 
##                             160                            3000 
## isUsingOneBigSeasonOnFewRecords 
##                               1
print(uStarTh)
##    aggregationMode seasonYear  season     uStar        5%       50%       95%
## 1           single         NA    <NA> 0.4541172 0.4188225 0.4564736 0.5286447
## 2             year       2020    <NA> 0.4801723 0.4488349 0.4876796 0.6412940
## 3             year       2021    <NA> 0.4280621 0.3888102 0.4258271 0.4512722
## 4           season       2020 2020001 0.4801723 0.3901247 0.4729320 0.6163569
## 5           season       2020 2020003 0.3993504 0.3877327 0.4484138 0.4839605
## 6           season       2020 2020006 0.2939407 0.2417530 0.3016328 0.3727611
## 7           season       2020 2020009 0.3376365 0.3284783 0.4120114 0.5594615
## 8           season       2021 2020012 0.4280621 0.3950345 0.4302720 0.4641630
## 9           season       2021 2021003 0.4280621 0.3742071 0.4258271 0.4512722
## 10          season       2021 2021006 0.3316664 0.2639192 0.3323809 0.3794716
## 11          season       2021 2021009 0.3632184 0.3288625 0.3656392 0.4183202
## 12          season       2021 2021012 0.4280621 0.3888102 0.4258271 0.4512722
# Define aggregation mode
EProc$sGetUstarScenarios() # by default, annual varying u* is used
##    season     uStar       U05       U50       U95
## 1 2020001 0.4801723 0.4488349 0.4876796 0.6412940
## 2 2020003 0.4801723 0.4488349 0.4876796 0.6412940
## 3 2020006 0.4801723 0.4488349 0.4876796 0.6412940
## 4 2020009 0.4801723 0.4488349 0.4876796 0.6412940
## 5 2020012 0.4280621 0.3888102 0.4258271 0.4512722
## 6 2021003 0.4280621 0.3888102 0.4258271 0.4512722
## 7 2021006 0.4280621 0.3888102 0.4258271 0.4512722
## 8 2021009 0.4280621 0.3888102 0.4258271 0.4512722
## 9 2021012 0.4280621 0.3888102 0.4258271 0.4512722
# Here, we use the single u* threshold
uStar <- uStarTh  %>%
  dplyr::filter( aggregationMode == "single") %>%
  dplyr::select( "uStar", "5%", "50%", "95%")
uStarDf <- cbind(season=na.omit(unique(uStarTh$season)), uStar)
EProc$sSetUstarScenarios(uStarDf)
EProc$sGetUstarScenarios()
##    season     uStar        5%       50%       95%
## 1 2020001 0.4541172 0.4188225 0.4564736 0.5286447
## 2 2020003 0.4541172 0.4188225 0.4564736 0.5286447
## 3 2020006 0.4541172 0.4188225 0.4564736 0.5286447
## 4 2020009 0.4541172 0.4188225 0.4564736 0.5286447
## 5 2020012 0.4541172 0.4188225 0.4564736 0.5286447
## 6 2021003 0.4541172 0.4188225 0.4564736 0.5286447
## 7 2021006 0.4541172 0.4188225 0.4564736 0.5286447
## 8 2021009 0.4541172 0.4188225 0.4564736 0.5286447
## 9 2021012 0.4541172 0.4188225 0.4564736 0.5286447
# if you want to use seasonal or annual varying u*:
# EProc$useSeaonsalUStarThresholds()# seasonal varying u*
# EProc$useAnnualUStarThresholds() # annual varying u*

# EProc$sPlotNEEVersusUStarForSeason(format = "pdf") # keep this line commented when use binder, then check the figure under /plots on github repo

Bonus training: create a figure comparing u* threshold determined using different aggregationMode.

06: MDS: gapfil NEE

“look-up” table (LUT)

Fluxes are expected similar if they are:

  • at similar environmental conditions

    • Rg \(\pm\) 50 \(W m^{-2}\), Tair \(\pm\) 2.5 \(°C\), and VPD \(\pm\) 5.0 \(h Pa\)
  • and close in time

    • increasing time window until enough observations

“mean diurnal course” (MDC)

Fluxes are expected similar if they are at the same time of the day and not too many days away;
The methods are effective for short gaps as the missing values are replaced by the average of response variables under similar weather conditions in a small-time window;

Quality flag increases with fewer variables and larger time windows

  • 0: true measurement
  • 1: gap-filled with good quality
  • \(>1\): gap-filled with lower quality
# fingerplot: inspect gaps needed to be filled
EProc$sPlotFingerprintY("NEE", Year = 2020); EProc$sPlotFingerprintY("NEE", Year = 2021)

# use MDS to gapfill flux data
EProc$sMDSGapFillUStarScens('NEE')
## Ustar filtering (u * Th_1 = 0.45411720625), marked 33% of the data as gap
## Initialized variable 'NEE' with 20978 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_uStar_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29068
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## .905
## Look up table with window size of 7 days with Rg
## ...................................................3092
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 20978, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.418822541913729), marked 31% of the data as gap
## Initialized variable 'NEE' with 20411 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_5%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29158
## Look up table with window size of 14 days with Rg VPD Tair
## ...........................................................843
## Look up table with window size of 7 days with Rg
## ..................................................3064
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 20411, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.456473581875), marked 33% of the data as gap
## Initialized variable 'NEE' with 21021 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_50%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................29062
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## .909
## Look up table with window size of 7 days with Rg
## ...................................................3094
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....529
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 21021, unfilled (long) gaps: 0.
## Ustar filtering (u * Th_1 = 0.528644729070933), marked 37% of the data as gap
## Initialized variable 'NEE' with 22024 real gaps for gap filling of all 35088 values (to estimate uncertainties).
## Full MDS algorithm for gap filling of 'NEE.Ustar_95%_fqc_0' with LUT(Rg, VPD, Tair) and MDC.
## Look up table with window size of 7 days with Rg VPD Tair
## ............................................................
## .............................................................
## .............................................................
## .............................................................
## .............................................................
## ...................................................28810
## Look up table with window size of 14 days with Rg VPD Tair
## ............................................................
## ...1076
## Look up table with window size of 7 days with Rg
## ....................................................3178
## Mean diurnal course with window size of 0 days: .
## ....................5
## Mean diurnal course with window size of 1 days: .
## ....................29
## Mean diurnal course with window size of 2 days: .
## ...................10
## Look up table with window size of 21 days with Rg VPD Tair
## ...................1
## Look up table with window size of 28 days with Rg VPD Tair
## ...................0
## Look up table with window size of 35 days with Rg VPD Tair
## ...................0
## Look up table with window size of 42 days with Rg VPD Tair
## ...................0
## Look up table with window size of 49 days with Rg VPD Tair
## ...................0
## Look up table with window size of 56 days with Rg VPD Tair
## ...................0
## Look up table with window size of 63 days with Rg VPD Tair
## ...................0
## Look up table with window size of 70 days with Rg VPD Tair
## ...................0
## Look up table with window size of 14 days with Rg
## ...................768
## Look up table with window size of 21 days with Rg
## ............681
## Look up table with window size of 28 days with Rg
## .....530
## Finished gap filling of 'NEE' in 27 seconds. Artificial gaps filled: 35088, real gaps filled: 22024, unfilled (long) gaps: 0.
# gaps in flux data are filled 
EProc$sPlotFingerprintY("NEE_50._f", Year = 2020); EProc$sPlotFingerprintY("NEE_50._f", Year = 2021)

07: Partitioning NEE into Reco and GPP

Reco is modeled and GPP is computed as Reco - NEE;

There are three approaches to partition FCO₂ available in REddyProc:

# TimeZoneHour: time zone offset from UTC (Coordinated Universal Time), without daytime saving
EProc$sSetLocationInfo(LatDeg = 46.2420, LongDeg =-89.3477, TimeZoneHour = -6)
EProc$sFillVPDFromDew() # fill longer gaps still present in VPD_f
EProc$sMRFluxPartitionUStarScens("NEE") # night-time
## Start flux partitioning for variable NEE with temperature Tair_f.
## Estimate of the temperature sensitivity E_0 from short term data: 84.51.
## Regression of reference temperature R_ref for 169 periods.
## Start flux partitioning for variable NEE with temperature Tair_f.
## Estimate of the temperature sensitivity E_0 from short term data: 84.51.
## Regression of reference temperature R_ref for 169 periods.
## Start flux partitioning for variable NEE with temperature Tair_f.
## Estimate of the temperature sensitivity E_0 from short term data: 84.51.
## Regression of reference temperature R_ref for 169 periods.
## Start flux partitioning for variable NEE with temperature Tair_f.
## Estimate of the temperature sensitivity E_0 from short term data: 84.51.
## Regression of reference temperature R_ref for 169 periods.
# If you want to expolore day-time or night-time partitioning
# EProc$sGLFluxPartitionUStarScens() # day-time
# EProc$sTKFluxPartitionUStarScens() # modified day-time
names(EProc$sTEMP)
##   [1] "sDateTime"         "Tair_orig"         "Tair_f"           
##   [4] "Tair_fqc"          "Tair_fall"         "Tair_fall_qc"     
##   [7] "Tair_fnum"         "Tair_fsd"          "Tair_fmeth"       
##  [10] "Tair_fwin"         "Rg_orig"           "Rg_f"             
##  [13] "Rg_fqc"            "Rg_fall"           "Rg_fall_qc"       
##  [16] "Rg_fnum"           "Rg_fsd"            "Rg_fmeth"         
##  [19] "Rg_fwin"           "VPD_orig"          "VPD_f"            
##  [22] "VPD_fqc"           "VPD_fall"          "VPD_fall_qc"      
##  [25] "VPD_fnum"          "VPD_fsd"           "VPD_fmeth"        
##  [28] "VPD_fwin"          "NEE_orig"          "NEE_f"            
##  [31] "NEE_fqc"           "NEE_fall"          "NEE_fall_qc"      
##  [34] "NEE_fnum"          "NEE_fsd"           "NEE_fmeth"        
##  [37] "NEE_fwin"          "season"            "Ustar_uStar_Thres"
##  [40] "Ustar_uStar_fqc"   "NEE_uStar_orig"    "NEE_uStar_f"      
##  [43] "NEE_uStar_fqc"     "NEE_uStar_fall"    "NEE_uStar_fall_qc"
##  [46] "NEE_uStar_fnum"    "NEE_uStar_fsd"     "NEE_uStar_fmeth"  
##  [49] "NEE_uStar_fwin"    "Ustar_5._Thres"    "Ustar_5._fqc"     
##  [52] "NEE_5._orig"       "NEE_5._f"          "NEE_5._fqc"       
##  [55] "NEE_5._fall"       "NEE_5._fall_qc"    "NEE_5._fnum"      
##  [58] "NEE_5._fsd"        "NEE_5._fmeth"      "NEE_5._fwin"      
##  [61] "Ustar_50._Thres"   "Ustar_50._fqc"     "NEE_50._orig"     
##  [64] "NEE_50._f"         "NEE_50._fqc"       "NEE_50._fall"     
##  [67] "NEE_50._fall_qc"   "NEE_50._fnum"      "NEE_50._fsd"      
##  [70] "NEE_50._fmeth"     "NEE_50._fwin"      "Ustar_95._Thres"  
##  [73] "Ustar_95._fqc"     "NEE_95%_orig"      "NEE_95%_f"        
##  [76] "NEE_95%_fqc"       "NEE_95%_fall"      "NEE_95%_fall_qc"  
##  [79] "NEE_95%_fnum"      "NEE_95%_fsd"       "NEE_95%_fmeth"    
##  [82] "NEE_95%_fwin"      "VPDfromDew"        "PotRad_5%"        
##  [85] "FP_NEEnight_5%"    "FP_Temp_5%"        "E_0_5%"           
##  [88] "R_ref_5%"          "Reco_5%"           "GPP_5%_f"         
##  [91] "GPP_5%_fqc"        "PotRad_50%"        "FP_NEEnight_50%"  
##  [94] "FP_Temp_50%"       "E_0_50%"           "R_ref_50%"        
##  [97] "Reco_50%"          "GPP_50%_f"         "GPP_50%_fqc"      
## [100] "PotRad_95%"        "FP_NEEnight_95%"   "FP_Temp_95%"      
## [103] "E_0_95%"           "R_ref_95%"         "Reco_95%"         
## [106] "GPP_95%_f"         "GPP_95%_fqc"       "PotRad_uStar"     
## [109] "FP_NEEnight_uStar" "FP_Temp_uStar"     "E_0_uStar"        
## [112] "R_ref_uStar"       "Reco_uStar"        "GPP_uStar_f"      
## [115] "GPP_uStar_fqc"

Bonus training: comparing different partitioning method;

08: Visualise the output

FilledEddyData <- EProc$sExportResults()
combined.df <- cbind(EddyData, FilledEddyData)
names(combined.df)
##   [1] "Year"              "Hour"              "Date"             
##   [4] "Month"             "DoY"               "NEE"              
##   [7] "LE"                "H"                 "Rg"               
##  [10] "Tair"              "Tsoil"             "RH"               
##  [13] "VPD"               "Ustar"             "TIMESTAMP_START"  
##  [16] "TIMESTAMP_END"     "PPFD"              "Tair_orig"        
##  [19] "Tair_f"            "Tair_fqc"          "Tair_fall"        
##  [22] "Tair_fall_qc"      "Tair_fnum"         "Tair_fsd"         
##  [25] "Tair_fmeth"        "Tair_fwin"         "Rg_orig"          
##  [28] "Rg_f"              "Rg_fqc"            "Rg_fall"          
##  [31] "Rg_fall_qc"        "Rg_fnum"           "Rg_fsd"           
##  [34] "Rg_fmeth"          "Rg_fwin"           "VPD_orig"         
##  [37] "VPD_f"             "VPD_fqc"           "VPD_fall"         
##  [40] "VPD_fall_qc"       "VPD_fnum"          "VPD_fsd"          
##  [43] "VPD_fmeth"         "VPD_fwin"          "NEE_orig"         
##  [46] "NEE_f"             "NEE_fqc"           "NEE_fall"         
##  [49] "NEE_fall_qc"       "NEE_fnum"          "NEE_fsd"          
##  [52] "NEE_fmeth"         "NEE_fwin"          "season"           
##  [55] "Ustar_uStar_Thres" "Ustar_uStar_fqc"   "NEE_uStar_orig"   
##  [58] "NEE_uStar_f"       "NEE_uStar_fqc"     "NEE_uStar_fall"   
##  [61] "NEE_uStar_fall_qc" "NEE_uStar_fnum"    "NEE_uStar_fsd"    
##  [64] "NEE_uStar_fmeth"   "NEE_uStar_fwin"    "Ustar_5._Thres"   
##  [67] "Ustar_5._fqc"      "NEE_5._orig"       "NEE_5._f"         
##  [70] "NEE_5._fqc"        "NEE_5._fall"       "NEE_5._fall_qc"   
##  [73] "NEE_5._fnum"       "NEE_5._fsd"        "NEE_5._fmeth"     
##  [76] "NEE_5._fwin"       "Ustar_50._Thres"   "Ustar_50._fqc"    
##  [79] "NEE_50._orig"      "NEE_50._f"         "NEE_50._fqc"      
##  [82] "NEE_50._fall"      "NEE_50._fall_qc"   "NEE_50._fnum"     
##  [85] "NEE_50._fsd"       "NEE_50._fmeth"     "NEE_50._fwin"     
##  [88] "Ustar_95._Thres"   "Ustar_95._fqc"     "NEE_95%_orig"     
##  [91] "NEE_95%_f"         "NEE_95%_fqc"       "NEE_95%_fall"     
##  [94] "NEE_95%_fall_qc"   "NEE_95%_fnum"      "NEE_95%_fsd"      
##  [97] "NEE_95%_fmeth"     "NEE_95%_fwin"      "VPDfromDew"       
## [100] "PotRad_5%"         "FP_NEEnight_5%"    "FP_Temp_5%"       
## [103] "E_0_5%"            "R_ref_5%"          "Reco_5%"          
## [106] "GPP_5%_f"          "GPP_5%_fqc"        "PotRad_50%"       
## [109] "FP_NEEnight_50%"   "FP_Temp_50%"       "E_0_50%"          
## [112] "R_ref_50%"         "Reco_50%"          "GPP_50%_f"        
## [115] "GPP_50%_fqc"       "PotRad_95%"        "FP_NEEnight_95%"  
## [118] "FP_Temp_95%"       "E_0_95%"           "R_ref_95%"        
## [121] "Reco_95%"          "GPP_95%_f"         "GPP_95%_fqc"      
## [124] "PotRad_uStar"      "FP_NEEnight_uStar" "FP_Temp_uStar"    
## [127] "E_0_uStar"         "R_ref_uStar"       "Reco_uStar"       
## [130] "GPP_uStar_f"       "GPP_uStar_fqc"
# Have a look at the gapfilled data
# diurnal pattern of FC
ggplot(combined.df, aes(x = Hour, y = NEE_50._f)) +
  geom_point(col = "grey") +  
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  stat_summary(fun = mean, geom = "line") + 
  facet_wrap(~ interaction(Month, Year)) +
  theme_minimal() +
  ggtitle("FC") +  ylab(expression(FCO[2]~'('*gC~m^{-2}~y^{-1}*')'))

# diurnal pattern of GPP
ggplot(combined.df, aes(x = Hour, y = `GPP_50%_f`)) +
  geom_point(col = "grey") +  
  stat_summary(fun = mean, geom = "line") + 
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +  
  facet_wrap(~ interaction(Month, Year)) +
  theme_minimal() +
  ggtitle("GPP") + ylab(expression(GPP~'('*gC~m^{-2}~y^{-1}*')'))
## Warning: Removed 14906 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 14906 rows containing missing values or values outside the scale range
## (`geom_point()`).

# diural pattern of Reco
ggplot(combined.df, aes(x = Hour, y = `Reco_50%`)) +
  geom_point(col = "grey") +  
  stat_summary(fun = mean, geom = "line") + 
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +  
  facet_wrap(~ interaction(Month, Year)) +
  theme_minimal() +
  ggtitle("Reco") + ylab(expression(Reco~'('*gC~m^{-2}~y^{-1}*')'))

Bonus training: gap-fill H or LE, and plot the diurnal patterns.

References